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Abstract We develop a Nitsche-based formulation for a general class of stabilized finite element 
methods for the Stokes problem posed on a pair of overlapping, non-matching meshes. By ex- 
tending the least-squares stabilization to the overlap region, we prove that the method is stable, 
consistent, and optimally convergent. To avoid an ill-conditioned linear algebra system, the scheme 
is augmented by a least-squares term measuring the discontinuity of the solution in the overlap 
region of the two meshes. As a consequence, we may prove an estimate for the condition number 
of the resulting stiffness matrix that is independent of the location of the interface. Finally, we 
present numerical examples in three spatial dimensions illustrating and confirming the theoretical 
results. 

Keywords Fictitious domain • Stokes problem • stabilized finite element methods • Nitsche's 
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1 Introduction 



Overlapping mesh methods offer many advantages over standard finite clement methods that re- 
quire the generation of a single conforming mesh resolving the full computational domain. With 
overlapping mesh methods, the computational domain may instead be described by a set of overlap- 
ping and non-matching meshes. In particular, different subdomains may be meshed independently 
and then collected to form the full domain. This feature is particularly useful in engineering appli- 
cations where meshes for physical components may be reused in different configurations. Another 
important example is the simulation of the flow around a complex object embedded in a channel. 
One may then create a mesh that discretizes a fixed and simple domain such as a cube or a sphere 
surrounding the complex object. This mesh may then be imposed on top of a fixed background 
mesh for the simulation of the flow around the object inserted at different locations in a domain 
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Fig. 1.1 The stabilized Nitsche overlapping mesh method presented in this work allows the simulation of Stokes flow 
around a complex object (here a simple two-dimensional airfoil) described by a matching mesh of its surroundings 
imposed on top of a non-matching fixed background mesh. 



of interest. A particular advantage of this approach is that it allows the creation of a fixed graded 
mesh to resolve boundary layers close to the surface of the complex object. This is illustrated in 
Figure [Ll] for a simple two-dimensional airfoil embedded in a channel. 

In this work, we introduce an overlapping mesh method for Stokes flow with constant viscosity 
across the artificial mesh interface. The Stokes problem reads: find the velocity u : fl C lR d — > M d 
and the pressure p : Q — ¥ M. such that 

-Au + Vp = f mf2, (1.1a) 
V • u = inQ, (1.1b) 
u = g on dn, (1- 1C ) 

where n denotes a bounded domain in R d , d = 2 or 3, with Lipschitz boundary dn, and where 



/ G L 2 (n) and g G H 1 ^ 2 (dn) are given functions. To satisfy (1.1b I, we assume that the mean 
value of g ■ n vanishes; n denoting the outward unit normal to dn. 

Our formulation is based on a general stabilized Galerkin finite element method for the Stokes 
problem and enforcement of the interface conditions via Nitsche's method. In order to prove 
stability, we let the least-squares stabilization terms extend to all elements that intersect the com- 
putational domain. As a result, these terms will appear twice in any overlap regions. In addition, 
we include a certain least-squares term penalizing the difference between the velocity solutions in 
the overlap region. This allows us to prove stability and optimal order error estimates as well as 
to control the condition number of the resulting algebraic problem. 

The method proposed here can be viewed as an extension to the Stokes problem of earlier 
work by Hansbo et al. |13j who developed a Nitsche overlapping mesh method for a second order 
elliptic model problem. Also, Becker et al. [5] presented a Nitsche extended finite element method 
for incompressible elasticity based on low-order ([Pf] d x Pq) elements. Moreover, the least-squares 
penalty of the velocity differences is related to the mesh tying approach proposed by Day and 
Bochev [7] , who formulate a least-squares problem for a system consisting of the partial differential 
equation together with the interface conditions. Note that in our method, the interface conditions 
arc enforced using Nitsche's method, while the least-squares terms on the overlap are only included 
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to prove the stability of the method and to control the condition number. In a related work |23j , 
we present a stabilized Nitsche fictitious domain method for the Stokes problem. 

The implementation of the overlapping mesh method in three space dimensions is a challenging 
problem. A realization of the method proposed in this work entails computing the intersection of 
arbitrarily superimposed tetrahedral meshes and integration over arbitrarily cut tetrahedra. Such 
a realization has been developed as part of the C++ library DOLFIN-OLM ( |http : //launchpad . | 
|net/dolfin-olm[ ) extending the FEniCS Project software [20] QH US El HU HI E] ■ For a discussion 
of the computational aspects, we refer to our previous work [22] and the related paper [23] . 

The remainder of this work is organized as follows. We first summarize our assumptions and 
notation in Section [2] The overlapping mesh method is then formulated in Section [3] Sections [4]|6] 
are devoted to the stability and a priori error analysis of the proposed method, while the condition 
number estimate is presented in Section [7] Finally, we demonstrate the proposed method for a 
sample application in Section [8] and present numerical convergence results and condition number 
estimates to support our theoretical results. 

2 Preliminaries 

In this section, we review the notation used throughout the remainder of this work. We also 
summarize a standard stabilized Stokes formulation to lay the foundations for the formulation of 
the overlapping mesh method in Section [3] 

2.1 Finite element spaces 

In what follows, H S (S7) denotes the standard Sobolev space of order s € N, defined on an open 
and bound domain fl with Lipschitz boundary dfl. We write (•, -) s _n^ II ■ \\s,n an d | ■ | s q for the 
inner product, norm and semi-norm on H s (fi), respectively. The index s will be dropped when 
s = 0. 

For a given, shape-regular tessellation T of f2, we let the associated discrete velocity space 
Vh C [H 1 (Sl)] d be the space of continuous piecewise polynomial IR^-valued vector fields of order 
k, and let the pressure space Qh C L 2 (Q) n C°(J7) consist of continuous piecewise polynomials of 
order I. To emphasis the order of the underlying polynomials, we occasionally write V£ and Q l h . 

2.2 Stabilized Stokes elements 

We recall the definition of consistently stabilized finite element methods for the Stokes problem, 
following Franca et al. [TU], Barth et al. [I]. We first define the bilinear and linear forms 

ah{uh,v h ) = {Vu h ,\7v h ) n , (2.1) 
bh(vh,Ph) = -(V • v h ,p h )a, (2.2) 
k(v h ) = (f,v h )n- (2.3) 

As categorized by Barth et al. [1], a wide class of consistently stabilized mixed finite element 
formulation for the Stokes problem can be recast in the following form: find (uh,Ph) € Vh, X Qh 
such that 

A h (u h ,p h ;v h ,q h ) = L h {v h ,q h ) \/{v h ,q h ) e V h x Q hl 

where 

A h (u h ,p h ;v hl q h ) = a h (u h ,v h ) + b h (v h ,p h ) + b h (u h ,q h ) - S h (u h ,p h ;v h ,q h ), 
L h (v, q) = l h {v) -5^2 h r(f, -otAv h + PVq h ) T , 

TET 
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and where the stabilization is given by 

S h (u h ,p h \v hl q h ) = S ^2 h\(-Au h + Vph, -aAv h + j3Vq h ) T - 

TET 

The parameters a and j3 are chosen from the sets {—1, 0, 1} and {—1, 1}, respectively and 6 denotes 
some positive constant. Barth et al. [4] point out that the choice (a,j3) = (1, 1) corresponds to the 
classical scheme introduced by Hughes et al. [16j . while on the other hand, the method by Douglas 
and Wang [5] can be constructed from the parameter choice (a,/3) = (—1, 1). In what follows, we 
will focus on these two families of stabilized methods for the Stokes problem. 



2.3 A domain decomposition model problem for the Stokes problem 



Let fl = (fix U Q2) be a domain in K d with Lipschitz boundary dfl, consisting of two (open and 
bounded) disjoint subdomains fix and fl 2 separated by the interface r = dflx Ddfl2- To develop a 
Nitsche based overlapping mesh method for the Stokes problem, we consider the following domain 
decomposition model problem for ( 1.1 ): find u : fl — > M 3 and p 



such that 



-Aui + Vpi = f i in Qi, i 

V • U{ = in i?^, i 

[u] =0 on r, 

[d n u — pn] =0 on r, 

u = on dQ. 



1,2, 
1,2, 



(2.4) 
(2.5) 
(2.6) 
(2.7) 
(2.8) 



Here and in the following, v, t — v\^. denotes the restriction of a function or vector field v to a 
subdomain fii. Furthermore, n is the unit normal of r directed from Q2 into J7i, d n u = n ■ Vu, 
and [v] = 1)2 — Vx denotes the jump in a function over the interface r. 




Fig. 2.1 Decomposition of the domain Q by introducing the artificial interface J\ The weak coupling along the 
interface _T by the Nitsche method allows independent meshes for the subdomains fix and Q2 , including overlapping 
meshes. 



The decomposition of J? into the two subdomains fix and fli motivates the introduction of the 
broken Sobolev spaces 



H s (f>xUfl 2 ) = H s (f2x)@H s (fl 2 ), 0, 



(2.9) 
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endowed with the norm 



The key idea in developing in Nitsche-type methods is to now replace the strong continuity 



conditions (2.6) and (2.7) by a weak formulation [231 HI 1121 EH 03] ■ This approach is analogous 
to that of discontinuous Galerkin methods for elliptic equations [3]. Starting with suitable finite 
element discretizations of H s (f2i) and 77 s (.f2 2 ), a weak formulation can be obtained by multiplying 
with test functions, integrating by parts and adding certain symmetrization and stabilization 
terms. A typical example of the resulting interface form will be given as part of the method we 
present in Section [3j 

We remark that the introduction of the interface T is purely artificial in our application case 
and solely serves the purpose of decomposing the domain into suitable subdomains to ease and 
decouple the subsequent meshing process. 



2.4 Overlapping meshes 



We consider a situation where a background mesh To is given for fl = U .f? 2 ) and another 



mesh 7~2 is given for the overlapping domain J? 2 (see Figure 2.1). Both meshes are assumed to 
consist of shape-regular simplices T. We note that the tessellation To of the background domain 
i? may be decomposed into three disjoint subsets: 

T = To,iUTo, 2 UT ,r, (2.10) 

where %,i = {T G To : T C T , 2 = {T G To : T C 7? 2 } and To,r = {T G To : \T H f2 t \ > 

0, i = 1,2} denote the sets of not, completely and partially overlapped elements relative to f? 2 , 
respectively. The meshes T\ and 71 are then defined by 

7T = T ,i U %,r, (2.11) 
T = {Tn7?i : Te 7?}. (2.12) 

Moreover, we introduce the tessellated domain [21 — [J Te -j-- T and the overlap region Q Q 



J7 2 n fi^. To ease the notation, we occasionally refer to Ti as 7^* and to f2 2 as J?|- ^ee Figure 2.2 
for an illustration of the various mesh parts and regions. Furthermore, we introduce the notation 
d e T* (the exterior facets) and d{TT , j = 1,2 (the interior facets), for the set of facets which 
belong to either one or two elements, respectively. Here, a facet means an edge in M 2 and face in 



M 3 . In accordance with (2.12), the sets d{T\ and d{T\ denote the corresponding set of intersected 
facets F D Q\ . 

In addition to the shape-regularity, we require that the mesh sizes are compatible over the 
interface T. More precisely, we assume that there exist a mesh-independent constant C > such 
that for all T G To and all T G 7i it holds 

C- x h T < hf < Ch T (2.13) 

whenever T n T D T ^ 0. Similar assumptions were made in Hansbo et al. [TB"] . 

Finally, we compose suitable finite element spaces for the overlapping meshes. For each mesh 
Tj* , £ = 1,2, let Vh,i be a finite element space of continuous fixed-order polynomials on T* ■ Similar 



to the broken Sobolev space (2.9), we introduce the composite finite element space for the whole 
domain ft 

Vfc=V M 0V fcl2j (2.14) 



and define the broken Sobolev space norm as before by 

\l,T 1 *UT 2 = \\ v h,l\\s,7? T II^MIIa.ra- 



v h\\l T *uT 2 = \\ v h,i\\i T * + hha\\l,T 2 - (2-15) 
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for Vh G W Note that both v^i and contribute to the norm in the overlap region Q . In 
particular, we define the composite finite element spaces Vh and Qh for the velocity and pressure, 
respectively, by 



Vh = V^ x v£ 2 , Q h = Q l h ,i QL,2 



(2.16) 
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3 A Nitsche overlapping mesh method for the Stokes problem 

Given a domain fl, overlapping meshes T\ and Ti and the composite finite element spaces Vh and 
Qh as introduced in (2.16), we define the bilinear form Ah by 



A h (u h ,p h ;v h ,q h ) = a h (u h ,v h ) + b h (v h ,Ph) + b h (u h ,q h ) 
+s h (u h ,v h ) - S h (u h ,p h :v h ,q h ), 



(3.1) 



where the bilinear forms ah, bh, Sh and Sh are given by 



a h (u h ,v h ) = (Vitft,Vw/ ( )fi 1 ux2 2 ~ ((d n Uh),[v h })r (3.2) 

- ({d n v h ),[uh])r + 7(/i~ 1 [it/ l ], [vh])r, 

bh(vh,qh) = -(V ■ v h ,qh)n l \jn 2 + (n • [vh], {qh))r, (3.3) 

Sh{u h , v h ) = (V(u h ,i - u hfi ), V(« h ,i - »fe,i))« , (3.4) 

S h {u h ,p h ;v h ,q h ) = 5 ^2 h 2 T (-Au h +^Jp h ,-aAv h + (3Vq h ) T . (3.5) 
Ter 1 *ur 2 

Here, (u) denotes a weighted average (v) = a±Vi + av2 with «i + ai — 1 for a field u defined on 
fl and discontinuous across the interface T. In accordance with Hansbo et al. |13j and to simplify 
the presentation, we choose (v) = v 2 - Finally, the linear form Lh is defined by 

L h (v,q) = (f,v)-S h 2 T (f,-aAvh + PVq h )T- (3.6) 

Ter 1 *ur 2 



With these definitions, the Nitsche based overlapping mesh method for the Stokes problem (1.1 1 
reads: find (uh,Ph) € Vh x Qh such that 

Ah{u h ,Ph\v h ,q h ) = L h (v h ,q h ) V(vh,qh) G Vfc x Qfe- (3-7) 



Remark 3.1 The forms s/,, and the mesh 7^* are defined such that contributions from both 
meshes are taken into account in the overlap region Q Q . These twice-counted contributions allows 
us to obtain stability and condition number estimates; this will be discussed further in Section [5] 
and [3 

We conclude this section by stating the following Galerkin orthogonality relation. 



Proposition 3.1 Let (u,p) £ [H 2 (Q)] d x H 1 ((l) be a (weak) solution of (1.1) and {uh,,Ph) £ 
Vh x Qh, the solution of the finite element formulation (3.7). Then 

A h (u-Uh,p-Ph;v h ,q h ) =0 V (vh,qh) 6 Vh x Q h . (3.8) 



Proof The proof is standard and follows from integration by parts, the continuity conditions ( |2.6[ ) 
- (|2~7|) and observing that s h (u, v h ) = = S h (u,p; v h ,q h ) for all (v h) q h ) € V h X Qh,. □ 



4 Approximation properties 

In this section, we establish the appropriate interpolation operators and provide interpolation 
estimates for use in Sections [5j|6] for the a priori error analysis of the method proposed in Section[3j 
With reference to Section [2~4} Vh will here denote some composite finite element space Vh = 
Vh,i © Vft.,2 comprising piecewise polynomial functions or vector fields defined on T* , i = 1,2. In 
particular, we shall let = or = Qh, in the subsequent sections. The constants C involved 
in the inequalities will depend only on fi ( , the regularity of the function spaces considered, and 
possibly the shape- regularity of Ti and the polynomial order of Vh,i, i — 1,2. 
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4.1 Norms 

Given the domain fl = Q\ U J?2 and corresponding overlapping meshes T* and T2, we introduce 
the following pairs of mesh-dependent norms for (v, q) £ [H 2 (f2)] d x H 1 (f2) and (v^, q^) £ Vh x 



(see (2.16) and (2.15)) 



= llVwl 



J?lUJ?2 



IVwIIt^ut, 

MlfliUfla' 
l?ll7?UT a i 



IKM 11-1/2,/^ 

KM ll-i/ w 



Ml 



ll/2,h,ri 
2 

1/2, /i,r' 



where 



and finally 



E 

Ter 2 



/ ' l T 2a |l' u llo,rnT! 



(f,g)IH2 



(4.1) 
(4.2) 
(4.3) 
(4.4) 

(4.5) 



(4.6) 
(4.7) 



Note that ||| • |||*-norms are defined on the regular meshes Tj* and T2 and therefore represent proper 
norms for the discrete finite element functions, while ||| • ||| is more suitable for functions defined 
on i? = Q\ U Q-x only. 



4.2 Trace inequalities and inverse estimates 



In the following, T* denotes always one of the regular, non-intersected meshes 7^* and 7^* . We 
first recall the following trace inequalities for v £ if 1 (J?) (or v £ [H 1 (£2)} d ) [13]: 



IMIar < C{h T 1/2 \\v\\ T + 4 /2 ||Vw|| T ) VT g T* 
Mlrnr < C(h~ 1/2 \\v\\ T + 4 /2 ||Vw|| T ) VT G T* 



(4.8) 
(4.9) 



The following inverse estimates for v £ V/, = V/ l (7"*) will also be frequently used. See [231 IT31 [T5] ) 
for proofs. 



h T \\Vv h \\ T ^ C\\v h \\ 
1/2 



V T e T* , 



^II^Hf <C||V^|| t VTeT*, 
/4 /2 ||d n ^||rnT < CHVwfcllr VT G T*. 



(4.10) 
(4.11) 
(4.12) 



We emphasize that hp in (4.12) (as in ( 4.11[ )) denotes the diameter of the entire facet F as part 
of the mesh T*. Finally, the estimation of the stabilization terms (3.5) will involve the inverse 
estimates 



J2 h 2 T \\Vv h f ^C\\v h f T , 



Ter* 



h 2 T \\Av h \\ 2 < C\\Vv h \ 



Vi'/, G Vh, 
Vt'/j G V h . 



(4.13) 
(4.14) 



Ter* 



To the end of this section, we state a Poincare type inequality for finite element functions on 
the overlapped mesh domain Ti*. To prove the proposition, we will need the following lemma, 
which is stated and proved in [2"5] : 
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Lemma 4.1 Let T = {T} be a mesh consisting of shape-regular elements T and take any two 
elements T\ and Ti sharing a face F. Furthermore, let u be a piecewise polynomial function defined 
on the macro-element M = T\ UT2. There exists a constant C depending only on the mesh quality 
parameters and the maximal polynomial order p of u restricted to each of the elements, such that 



\H 2 Tl *s c(H|| 2 + ]T4 3+1 ([d>], [diu]) F ). 



(4.15) 



Also note that the inverse inequalities (|4.10|) and (|4.1l|) can be generalized, by the standard scaling 

■u\\t for i j, 



argument, to inverse estimates of the type 



h{\\D ] u\\ T Ch^WD^-' 1 " 
h^Wd'iuWp^CWD^uWT, 

where D^u denotes the j-th total derivative of u. 



(4.16) 
(4.17) 



Proposition 4.1 Let fl\ be the overlapped domain and T\ be the overlapped mesh, and let Vh,i 
be a finite element space on 7^* consisting of continuous piecewise polynomials. Then 



IKIIt? < C{\\v h \\ 2 ni + Y, h 2 T \\Vv h \\%) 
\\v h f T * < C(K|| 2 fii + ||V^,), 



for all v h £ Vh,i- 
Proof By definition 



\Mh - E 

TeTo- 



\vhWr 



E 



KIItv 



Applying Lemma 4.1 and subsequently invoking ( |4.17 1 and (4.16) give 



E < C E (IKIIt + h%\\Vv h 

TiGTo.r TETf 



(4.18) 
(4.19) 

(4.20) 
(4.21) 



where the shape-regularity of 7^* provides a bound on hp in terms of hx- Combining (4.201 

and ( |4.21| yields ( |4T8| |. 

Since hp ^ l^i I, the second Poincare type inequality (4.19) is a simple consequence of the first 
one. □ 



4.3 Interpolation estimates 

The aim of this section is to construct and analyze an interpolation operator 7T/j : L 2 (fi) — > Vh and 
a corrected interpolation operator 7r£ : [H 1 (Q)] d — > Vh- The construction is based on extending 
the standard Scott-Zhang interpolation operator |26j . 

First, for s ^ 0, there exists a linear extension operator £ s : H s (Sli) —> H s (f2l) and a constant 
C > such that £ s v\q 1 = v and 

\\£'v\\. t at < C\\v\\ Si0l (4.22) 

for all v G H s (nx) [37]. 

Next, for i = 1,2, let 717^ ■ L 2 (S1*) —> Vh.i be the Scott-Zhang interpolation operator and 
recall the standard interpolation error estimates |26j : 

Wv—tt^v^t ^Ch s - r \ Vl \ sMT)l 0<r<s VT e 77, (4.23) 

\\v i -n h , i v i \\ r , F ^Ch s - r - 1 / 2 \v i \ sMF) , 0<r<s VFedtT*, (4.24) 

for all vi £ H s (£l*) for i = 1,2. Here, w(T) is the patch of all elements sharing a vertex with T; 
the patch oj(F) of a face F is defined analogously. 
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Altogether, we define 



s > 1, 



by 



(4.25) 



(4.26) 



The estimate (4.22) together with the interpolation error estimate (4.23) for the Scott-Zhang 



interpolation operator imply the following interpolation estimates: 

\\v - 7r h v\\ rt T ^ Ch s - r \v\ sMT) , O^r^s VTeTiUTj, (4.27) 

\\v - ir h v\\ r , F < Ch s - r -^ 2 \v\ sMT} , O^r^s VF G d{T\ U dffi- (4.28) 

We now return to our specific finite elements spaces Vh and Qh- We continue writing 7T/j for 



both the interpolation operator ( |4.25[ ) defined for [H k+1 (f2)] d Vfi and H l+1 (f2) Q l h . The 
following lemma provides an interpolation error estimate in the ||| • ||| norms. 

Lemma 4.2 There is a constant C > such that for all v e [H k+1 (Q)] d , q e H l+1 (Q) it holds 
that 



\\\v - n h v\\\ ^ Ch k \v\ k+1> a, 
\\\q - ir h q\\\ ^ Ch l+1 \q\ l+1>Q , 
(v - n h v, q - TT h q)\\\ ^ C{h k \v\ k+1 . n + h l+1 \q\ l+1 j 2 ). 



(4.29) 
(4.30) 
(4.31) 



definition of n^, the interpolation estimate (4.23) and the continuity of the extension operator 



Proof For the proof of (|4.29|) we refer to Hansbo et al. [To]. The estimate (|4.30|) follows from the 

or 
□ 



We conclude the section by referring to a modified Scott-Zhang interpolation operator intro- 
duced by Becker et al. [5] to establish a Nitsche extended finite element method for incompressible 
elasticity. There, the standard Scott-Zhang interpolant was corrected to gain some additional or- 
thogonality properties on the interface. Although only used for continuous piecewise linear velocity 
fields in combination with piecewise constant pressures in Becker et al. [5], we will show in the 
next section that the interpolant can be utilized to prove a sufficient inequality for bh{-, •) adapted 
to the case of overlapping meshes. 

To construct the corrected interpolant, Becker et al. [5] regrouped the elements T G 7o,_r into 
a collection r u of patches u>i in such a way that each patch can be associated with a hat function 
<I>i constructed from the finite element functions defined on elements Tcwj. Moreover, taj and <Pi 
have the following three properties: there are constants c, C > such that 

(i) ch diam(wi) ^ Ch 

(ii) ch < /rn<* *< < Ch 

(hi) ch- 1 < |V<?i(x)| < Ch- 1 

Loosely speaking, these properties guarantee that patches and hat functions behave as standard 
elements and shape functions. Although Becker et al. [5] only elaborate on the construction in 
two dimension and for a non-fitting interface, a similar construction is possible in any spatial 
dimension d and for the case where the interface is defined by the boundary of another mesh. 

By adjusting the degrees of freedom associated with the hat functions an interpolation 

operator 

can be defined that satisfies an orthogonality property of the form 



and an interpolation estimate of the form 

IK -Afc.i'Willo.flj + Ml V(vi - 7r, c l! i'Wi)||o,r2* «S C7i||Vv||o,rt*. 



(4.32) 



(4.33) 
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Finally, we can construct a corrected Scott-Zhang interpolation operator 7r? for the composite 
velocity spaces Vh by 

«MffW-» VuQVw (4 34) 

7T h( v ) = "m^i ®ir h ^v 2 

where Tth,2 is the standard Scott-Zhang operator on f2 2 - 

We will need the following continuity property of the corrected interpolation operator 7r£ with 
respect to different norms. 

Lemma 4.3 Let v E There exists a constant C > such that 

|||7r£v|||* < C||«||i,o lU fl a . (4-35) 

Proof By definition, 

IK^III* = \\^> h \\^, uT2 + \\{d n K c h v h )\\ 2 _ 1/2hr + \\Kv h }\\ 2 1/2hr . 

The estimate for the first term follows directly from the definition (and boundedness) of 7r£ and the 
continuity of the extension operator 8. Using the inverse estimate (4.12), the bound for the second 
follows from the bound of the first. To estimate the last term, we first note that ||[7r£i>]||i/2.h,r = 
|| [tt^v — v]\\i/ 2 .h,r since [v] = for v G i7 1 (J7) and thus 

\\K v ]\\ 2 i/2,h,r ^ \K V - v i\\l/2,h,r + \\< v - v 2\\l/2,h,r- 



Another application of the trace inequality (4.9) combined with interpolation estimates (4.33) 
yields for i = 1,2: 



\rnT 



Wh v ~ v if-i/2,h,r < h T 1 hh v - v \ 

< ]T h-^h-^v-vW^ + hWV^v-v)^) 

TGT* 

< \\v\\l T . < C\\v\\l niUn2 , 



□ 



5 Stability estimates 

In this section, we first prove that the form ah + Sh is continuous and coercive with respect to the 



norm in Proposition 5.1 that by t is continuous with respect to both the 



and the 



norms in Proposition |5.2| and an additional inequality for bh in Proposition |5.3| Subsequently, these 
estimates are used to show that the complete stabilized Nitsche form Ah satisfies the Babuska- 



Brczzi inf-sup condition in Theorem 5.1 



We begin by demonstrating that ah + Sh is continuous and stable with respect to the norm 



Proposition 5.1 Let ah and Sh be defined by (3.2) and (3.4). There are constants c > and 
C > such that 



a h (uh,v h ) + s h (u h ,v h ) < C\\\u h 

|2 



\Vh\\\* Vu h ,v h e V h , (5.1) 

< a h (v h , v h ) + s h {v h ,v h ) Vv h eV h . (5.2) 



Proof A proof in the absence of the term Sh and with 



Vu h ,v h e V h , 
Vv h e V h . 

replaced by ||| • ||| is given in Hansbo 



iVt^H^ the estimates (5.1) and (5.2) follow 



et al. [15] . Since by definition |||t)/,|||, = 

directly from the equivalence || Vvh\\j-*uT2 ~ ll^^llfiiu^ + ll^ 7 (' u 'i.i — ^ /i , 2 ) 1 1 ?2 and the definition 
ofs h . □ 
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Next, we state the continuity properties of bh with respect to the various norms. 



Proposition 5.2 Let bh be defined by (3.3). There exists a constant C > such that 

V(v,q) G [i? 1 (X?)] 3 x L 2 (f2), (5.3) 

(5.4) 



b h (v,q)^C\\\v\\\\\\q\\\ 
bh{v h ,qh) < V(vh,<lh) €V h x Q h . 



Proof The first inequality is a trivial consequence of the fact that the interface term [v] vanishes 
for v G [i/ 1 (i7)] 3 . To prove the second estimate, it is enough to estimate the interface term 
(n ■ [vh], (qh))r- Recall that we chose (qh) = qn,2- We combine the Cauchy-Schwarz inequality 
with p~9] ), ( |4T0[ ) and ( |4.12[ ) to obtain 

((q h ),n- [v h ]) r ^ \(qh,2,n- v ht i) r \ + \(qh,2,n- v h ,2)r\ 

^ lkft,2ll-i/2,/i,r {\\vh,i\\i/2,h.r + \\vh,2\\i/2,h,r) 
<C\\q h!2 \\r 2 (\\Vv hA \\ T . + \\Vv h!2 \\ T2 ) 

< CGIk/tllWII^III*)- 

□ 

The proof of the inf-sup condition for stabilized Stokes elements often involves an inequality 
known as Verfiirths trick [28]. In the next proposition, we present a version adapted to the case 
of overlapping meshes. 

Proposition 5.3 There are constants C\ > and C 2 > such that for all qh G Qh: 

sup irr^^w^-^t J2 4liv^|||) 1/2 . (5.5) 

v„ev h \\\Vh\\\* TeT 1 -ur 2 

Proof Let qh G Qh be given. The div-operator maps [H 1 (f2)] d onto L 2 (fl) and there is a constant 
C > such that ^ C|| divw||j7 [TT] for all v. Choosing v G [i? 1 (J7)] rf such that divu = —qh 

thus gives 

bh(v,qh) = \\qh\\l^C- 1 \\qh\\n\\v\\ 1 ,a, (5.6) 

where we have also used the fact that [v] = 0. 

Next, we take Vh — ^f t v, where 7r£ is the corrected Scott-Zhang intcrpolant introduced in 
Section l4~3l It follows that 

bh{v h ,qh) = b h {-K c h v,qh) = b h (TT c h v - v, q h ) + b(v, q h ) (5.7) 

^ b h (ir c h v -v,q h ) + C- 1 ||<2,J r2 |M|i, r 2- (5.8) 

To proceed further, we recall that (qh) = qh,2, n = n 2 and observe that 

2 

' i^h^i - Vi),qh,i)r = (n • [tt,> - v] 1 q h , 2 )r + [n ■ (n^Vx - v{), [q h ])r- 



Using this identity and integrating the first term in (5.8) by parts, we derive that 
h(^h v - v >1h) = -(V • (TT c h v - v),q h )n lU n 2 + (n • [n c h v - v],q h , 2 ) r 



i=l TeTi 



ii 



Rewriting (iriv — v,\7q h 



Tnf2 



(/i T 1 (7r£i; — v),hTyqh)mn and applying the interpolation esti- 



mates (4.27) and (4.33), the first term / can be bounded from below by 



\I\>-\\Vv\\ T * UT2 { ]T 4|| Vg, 

TGT 1 *UT 2 



,2\l/2 
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To estimate the second term //, we exploit the orthogonality property (4.32) and write 



: .r u 



(5.9) 



where qh,i\u = t^t f u qh,i da;, i = 1,2. Now each contribution in the sum (5.9) can bounded by 



(n • (7r£ 4 t>i - gv - qh,i)rnu 

^ ll n ' ( 7r h,i v i ~ v i)\\i/%h,rnu>\\qh,i ~ <ih,i\\-i/2,h,rn« 



(5.10) 



Using again the trace inequality (4.9), the last factor can estimated by 

\\qh,i ~ qh,i\\-i/2,h,rnu ^ C\uj\\\ Vqh,i\\u < Ch T \\Vqh,i\\u, 

since g/^j — g^j has mean value zero on u and so \\qh,i — <jh,i\\ui ^ M HVg/^jHw. To arrive at 
an estimate for first factor in (5.10), we use the trace inequality (4.9) in combination with the 



interpolation estimate (4.33) and the continuity of the extension operator: 



Thus 



\\ n ■ (K.,i v i - ui)||i/2,ft,rnu) < llVwil^nfii- 



(n • - «i), [g/J>| < 51 MI|Vg/alU + W^qhMW^lW^fh 

<c( #livrf) 1/9 ||vt7|| ni . 

Ter 1 *uT 2 



Putting / and II together, we see that 



\b h ffiv-v,q h )\^C{ 4l|V%|||) 1/2 ||V«||n, 

TgT 1 *uT 2 



which together with ( |5.8| ) gives the estimate 



v U2 



/2 



TGT 1 *UT 2 



Now the estimate follows from the continuity (4.35) of the operator 7r 



Remark 5.1 Although looking very similar to the standard "bad- inequality" , the inequality (5.5) 



bears some subtle differences. Note that we count contributions HVg^llr twice in the overlap 
domain Qq- This is mainly a consequence of the trace inequality (4.9). The L 2 -norm of qh on the 



other hand is only taken on fl\ U fl^\ thus requiring special consideration in the subsequent proof 
of the inf-sup estimate for A% using the norm |||<7/i|||* = |] Q/i || 7^*uTi ■ 

We conclude the section by proving the inf-sup condition for the bilinear form Ah with respect 
to the norm |||(-, 



Theorem 5.1 Let Ah be defined by (3.1 )-(3.5 ) and assume that either {a, /3} = {1, 1} or {a, /?} = 
{— 1, 1}. Then there is a constant c such that 



Ah(u h ,Ph]Vh,qh) . ,,,, s ln u/ , , T/ „ 

sup — 77T- w . ^ c\\\{u h ,Ph)\\\* V(u h ,Ph) eV h xQ h . 

(v h ,q h )eV h xQ h \\\(Vh,qh)\\\* 



(5.11) 
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Proof The proof follows the presentation given in Franca et al. [TO]- Let (uh,Ph) be given and 
consider the case {a,/?} = {1, 1}. 

Step 1: Choose = —w such that the suprcmum in (5.5) is attained. By scaling w, we 
can assume that |||w|||* = ||p/t||r2 1 ur2 2 - Inserting {vh,qh) = (— w,0) into Ah, we obtain via (5.1| 



and (5.5 



A h (u h ,p h ;v h ,q h ) = -a h (u h ,w) - s h (u h ,w) - b h (w,p h ) + S h (u h ,p h ;w,0) 
^ —C a \\\uh\\\*\\ph\\n 1 un 2 
+ {Cx\\p h \\n^n 2 -C 2 { ]T h T\\Vq h \\ 2 T ) 1/2 }\\ Ph \\ 

TeT 1 *uT 2 

+ 6 h 2 T (-Au h + \7 Ph ,-Aw) T ■ (5.12) 

Ter 1 *ur 2 



1/2 



We estimate / using Cauchy-Schwarz and the inverse inequality (4.14|: 

/ > -Cj5\\Vw\\ t ^t. (l|Vw,|| Tl *ur 2 + ( ]T 4I|Vp,Jt) 

Ter 1 *ur 2 

>-Cj*||p h ||« lU « a (lKIIU + ( E 4llv^|| T ) 1/2 ) 

TeT 1 *UT 2 

where the last inequality above follows as || Vo||7-*u7i ^ IIMII* f° r all v an d |||' u '|||* = ||p/i||fiiUf2 2 - 
Inserting this final estimate for / into (5.12| and using the inequality ab ^ ea 2 + (4e) _1 6 2 for any 
e > 0, we arrive at 

A h (u h ,p h ; -io,0) > (Ci-e(C a + C 2 + 5C I ))\\p h \\j hun2 
-(C a + ^ / )(4 £ )- 1 ||K||| 2 

-{c 2 + 8c I ){^r i h T\\vph\\ 2 T 

T&T 1 *uT 2 

>c 1 \\p h \\ 2 niuna -c 2 \\\u h \\\l-c 3 £ 4liv^||| 

TeT 1 *uT 2 

with positive constants Ci, C2, C3 for an appropriate choice of e. 

Step 2: Now choose (vh,qh) = {uh,—Ph}- Proposition 5.1 and subsequently (4.14) give that 

A h (u h ,p h ;u h ,~p h ) = a h (u h ,u h ) + s h (u h ,u h ) - S h (u h ,p h ;u h , -p h ) 

> c a lll-ix^lH^ - 5 X h^(-Au h + Vp h ,-Au h -Vp h ) T 

rer 1 *ur 2 

>(c„-«C , i)|||« A |||» + « X ^tII^IIt 

Ter 1 *ur 2 

TGT/UT2 

with positive (7 as long as < S < c a Cj l . 

Step 3: To complete the proof, we combine step 1 and step 2 by taking (vh,qh) = (uh—rjw, —ph) 
for some 77 > 0. By choosing rj sufficiently small, 

A(uh,Ph]Vh,qh)>Ci\\\u h \\\l + C 2 \\ph\\a lU a 3 +Ca X ^IIV^Hr 

TeTi*uT 2 

^C 1 ||K|||2 + 5 2 ||K|||2> m m((7 1 ,C 2 )|||(i ife ,^)|||2, 
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for some constants C\ , C 2 , C3 > where we used Proposition |4.1| to find that for some constant 
C 2 > 

C2\\Pk\\l lU n 2 +C 3 h T\\VPh\\ 2 T > C 2 \\\ Ph \\\l 

TeT 1 *UT 2 

+ »7|||u>|||* ^ (1 - 



(u h ,p h )\\\*, the estimate (5.11) follows 



Since \\\(v h , q h )\\\* < \\\(u h ,p h )\ 

The proof for {a, (3} = { — 1, 1} differs only in the derivation of the final estimate in step 2 
and can be taken from Franca et al. [10] with the presented adaption to the overlapping meshes 
formulation. □ 



Remark 5.2 We would like to comment on the role of the different ghost-penalties. As remarked 
earlier, the bad-inequality 5.5 "mixes domains" in the sense that it the contributions ||V(7/j||t 



are taken from the elements in 7f* U T 2 , so that the contributions from the overlap region f2 Q 
are counted twice, while \\qh\\ is evaluated only once on fl\ U fl 2 . The role of the pressure terms 
in the least-squares ghost-penalty X)t<eT*uT 2 ^ri~^ u h + Vp/,, —Av^ + Vqh)r is two-fold. First, 
they compensate the negative contributions — X^tst* ^tII^^IIt m the bad- inequality. Secondly, 
they allow in combination with Lemma 4.1 to pass from Hg/iH^ur^ 



to pass from \\q h \\n 1 yjn 2 to \\qh\\T*uT 2 - We further 
note that the velocity terms in the least-squares ghost-penalty make it necessary (via the inverse 
inequality ( 4.14 1) to control HVw^il 
Uh,2), V(f M - v h . 2 )) n 



\a o . This is precisely the role of s/ l (it^,t)/ l ) 



(V(«A, 



6 A priori error estimate 

The previous results on the interpolation errors, the Galerkin orthogonality of the discretization 
and its stability enable us to state the following a priori estimate for the error in the discrete 
solution. 

Theorem 6.1 Let k ,1 ^ 1 and (a, (3) = (±1, 1). Assume that (u,p) € [H k+1 (Q)] d x H l+1 (Q) is a 
(weak) solution of the Stokes problem (1.1 ). Then the finite element solution (xih,Ph) € Vu x Ql 
defined in (3.7 1 satisfies the following error estimate: 



\\\(u - u h ,p-p h )\\\ ^ C (h k \u\ W2 + h l+1 \p\ l+li a) ■ 
Proof By applying the triangle inequality, we have 

|||(tt - u h ,p-Ph)\\\ < Ml (1* - TT h u,p- TThP)\\\ + WK^hU - u h ,ir h p-p h ) 



(6.1) 



Given the interpolation estimate (4.31), it is suffices to bound the discretization error (71^ « 
UhiKhP — p)- By Theorem 5.1 there exists a (vh,Ph) such that |||(»fc, jJ/»)|||* ^ 1 and 

c\\\(ir h u-Uh,Trhp-ph)\\\* < Ah(uh - n h u,p ~ n h p;v h ,q h ). 



Thus, the Galerkin orthogonality (3.8) and the definition of Ah give 

c| I - u h ,7r h p - Ph)\\\* ^ A h (u- nhU,p- ir h p; v h ,q h ) 

= a h (u - ir h u,v h ) + s h (u - Ti h u,v h ) 



+ b h (v h ,p- n h p) + b h (u - TT h u, q h ) 

II /// 
-S h (u - ir h u;p- ir h p,v h ,q h ) . 



IV 



To estimate the first term, we use (5.1) and the interpolation estimate (4.29) to obtain 
|/| < C|||ti-7r h w|||*|||»fe|||» < Ch k \u\ k+1 ,«||K|||* < Ch k \u\ k+ha . 



16 



Andre Massing et al. 



Recalling that (p) = p2 , the second term 77 can be estimated by 

|77| = |(V ■ v,p- ■K h p)n l un 2 + (n ■ [«], (p - ^hP))r\ 

< llblll IIIp-^pIII + IIMIIi/2,/i,r||P2 - 7Th,2P2||-i/ 2 ,h,r 

< llblll IIIp-^pIII +C||[«]||l/2Ar(l|P2 -7Tft,2P2||x2 2 + ^11 V(p 2 - 7T ft , 2 P2) II f2 2 ) 

^C\\\v\\\h l+1 \p\ l+1 j 2 ^Ch l+1 \p\ l+hn , 



where we used the trace inequality (4.8) in the penultimate step and finally (4.10) in combination 
with the interpolation estimate (4.30). The third term can be treated similarly to obtain 

\iii\ < ch k+1 \u\ k+1>n \\\q h \\u < ch k+1 \ u \ k+U2 . 

For the last term IV, several applications of Cauchy-Schwarz yield 

\IV\=5\ h 2 T {-A{u-it h u) 1 +V{p-n h p) 1 ±Av h +Vq h ) T \ 
Ter 1 *ur 2 

<C( '4ll^(«-^«)ll7?ur» + ^'ll v (P- 7r ^)llW a ) 1/a 



Ter 1 *ur 2 

X( h T\\Av h \\ 2 T , ur2 +h 2 T \\Vq h \\ 2 T » uT2 ] 
TeT 1 *UT 2 



1/2 



After applying the inverse inequalities (4.14) and (4.13) to the second term, the interpolation 
estimate (4.27) implies 

\IV\ sc C{h k \u\ k+1J2 + fc ,+1 |p|i + i,fl)|||(t; h) flf h )||| < C(h k \u\ k+ i, n + h l+1 \p\ l+hn ). 

Summing up, we obtain that 

|||(u - u h ,p- Ph )\\\ sc |||(« - n h u,p- 7r hP )\\\ + \I\ + \II\ + \III\ + \IV\ 



^C(h k \u\ k+1 , n + h l+1 \p\ l<n ), 



which concludes the proof. 



□ 



7 Condition number estimate 

To conclude the analysis of the discretization presented in Section [3j we here show that the 
condition number of the associated stiffness matrix is uniformly bounded by Chr 2 independently 
of the location of the overlapping mesh 7a • The proof of the condition number estimate follows 
the approach of Ern and Gucrmond [9] . 

Let (v h , q h ) = Yli=i Vifi where {tfi\f=i is a basis for the element space 14 x Q h ; V = {Vi}i thus 
denotes the expansion coefficients of (vh, qh) in the given basis. Similarly, let W label the expansion 
coefficients of fields denoted (w h ,r h ). Denote the inner product in M. N by (V, W)n = J2iLi ViWi 
and the corresponding norm by |V|jy = (V,V)n- The stiffness matrix A associated with the 



form (3.1 1 is then defined as: 

(AV, W)n = A h (v h ,q h ;w h ,rh) Vv h ,w h e V h ,Vq h ,r h e Q h . (7.1) 

Since essential boundary conditions are applied for the velocity on the whole of d£2, the discrete 
pressure is only determined up to a constant, and so the matrix A has one zero eigenvalue. 
Throughout the remainder of this section, we therefore instead interpret A as A : R N / ker(A) — > 
ker(A)' L , which is a bijective linear mapping by definition. The condition number of the matrix A 
is then defined by 

k(A) = I^M^Iiv, (7-2) 
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IV 



with 



\A\ 



N 



\Ax\ N 
sup -r-j — • 

igR«\0 \ X \ N 



(7.3) 



We now state some estimates that will be needed in the derivation of the estimate of the 
condition number. First, for a conforming, quasi- uniform mesh T with mesh size h and a finite 
element space Vh defined on T, it is well known that there are positive constants cm and Cm 
only depending on the quasi-uniformness parameters and the polynomial order of Vh such that 
the following equivalence holds: 



cm 



h d/2 \V\ N < \\v h \\ < C M h d ' 2 \V\ N Vv h G V h - 



(7.4) 



A second ingredient is an inverse estimate and a Poincare inequality in the appropriate norms 
which we state in the following two lemmas. 



Lemma 7.1 There is a constant Ci such that 

\\\vh\\\* «S C I h- 1 \\v h \\ T * uT2 , 



K,gOIII*<^~lK,<Zfc)lk*u 



Vuft G V h , 

V(v h ,q h ) eV h x Q h . 



(7.5) 
(7.6) 



Proof By definition 



\v h \\\l = HVcfcH^uTj + ||<^t>fc)||l 1/aAr + ||[«fc]||? /aA r- 



Applying the standard inverse inequality (4.101 on each mesh T\ and 72 separately, the first term 
is bounded by Ch~ 2 \\vf l \\j-« u7 - 2 . A similar bound can be derived for the second term by (4.111 
and (4.10). The estimate for the remaining term follows in the similar manner by a combination 
of (|4.9[) and again ( |4.10| ). 

The second estimate (7.6) is an immediate consequence of the first and the fact that 1 = C ^ 
h^ 1 diam(J?). □ 

Lemma 7.2 There is a constant C such that 

\\v\\t?ut 2 < CpIIHH*, 
Proof By Proposition |4.1| we have 

Hl 2 ) 



Vv G V h . 



Mtcut, = IMIrr + IMIr 2 < + C(\\Vv\\* T . 



(7.7) 



<<7(|||«|||? 



I J?iUJ?2 



), 



for all v G Vh- Hence, it remains to derive a bound for a llfll^ui^- 

We use a duality argument to estimate ||i>||r2iUJ? 2 - Let <p G H 2 (f2) be the solution of the dual 
problem —A<p — v with boundary conditions <p — on df2. Multiplying with Vh and integrating 
by parts we get by using the Cauchy-Schwarz inequality and the trace inequality (4.9) for <fi: 

TnOi 



I r?iUJ? 2 



TeTi 



- ^2 (Vf, V0) T - (M,n ■ V0)r 

TeT 2 




V«l 



/i 2 |V0|? T ) 



Finally, using standard elliptic regularity 6 ||</>||2.r2 ^ \\ v \\n and division by ||i>||r2, the estimate 
follows. □ 
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Theorem 7.1 The condition number of the stiffness matrix satisfies the estimate 

k(A) ^ Chr 2 . (7.8) 



Proof Recalling the definition of the condition number in (7.2), the proof consists of deriving 

ve use tl 

(V, W) N 



estimates of \A\n and \A 1 |at. In the following, we use the following well-known equivalence for 
the Euclidean norm \x\n- 



\V\ N = sup 

w^o \W\n 



Estimate of \A\n- We have by the definition of A 

(AV, W) N 

\W\ N 
A h (v h ,q h ;w h ,r h ) 



\AV\ N =sup 

W^O \W\N 



= sup 
< C A sup ■ 



\w\ N 



{w h ,r h )\\L 



< CaCj sup 

< C A CjC 2 M \V 



\w\ N 

h~ l \\(v h ,q h )\\ T 'yjT 2 ' h- 1 \\(w h ,r h )\\ r * uT2 



\W\ 



N 



N • 



(7.9) 

(7.10) 

(7.11) 

(7.12) 
(7.13) 



Here, inequality (7.11 ) follows by the continuity of Ah with respect to the norm |||- 1||*, (7.12 ) follows 



|V|jv and taking the supremum over all V ^ 0, we conclude that 



by using the inverse estimate (|7.6|) twice, and finally (|7.13|) results by applying (7.4). Dividing by 

(7.W) 



\A\ N < C A CfC 2 



Estimate of \A~ 1 \n- Starting with (7.4) and sequentially using the Poincare inequality ( |7.7[ ) (ex- 
tended to the product space), the inf-sup stability (5.11 ) of Ah and finally the Poincare inequality 
again, we arrive at the following chain of estimates: 

h d/2 \V\ N < c^\\(v h ,q h )\\ T *uT, *S c£C P \\\(v h ,q h )\\\. 

^°M^ C A ||| (wk , rfc) |||, - C M^PC A ||| (wfe)r/i) || U 

^c M C P c A \AV\ N h mWhM 
^c^GlcA 1 \AV\Nh- d ' 2 

Setting V = A~ X W we have W = AV and the inequality now reads l^" 1 !^^ < C 2 P c^ 1 h ~ d \W\ N 
for all W ^ 0, or in other words 



\A- 1 \ N K t c^ti i P c^hr d . 



Finally, combining (7.14| and (7.15) yields (7.8). 



(7.15) 
□ 
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8 Numerical examples 

We conclude this paper with three numerical tests; all in three spatial dimensions. The numer- 
ical experiments were carried out using the DOLFIN-OLM library (http://launchpad.net/ 
dolf in-olm). We first corroborate the theoretical a priori error estimate with a convergence 
experiment. Second, we numerically examine how the location of the overlapping mesh in relation 
to the background mesh affects the condition number. Finally, we demonstrate how the method 
presented and the features provided by DOLFIN-OLM can be applied to the flow around an airfoil 
in a three dimensional channel. The experiments were performed with Vh = V^, Qh = Q\ and 
,3=1. 



.1 Convergence test 



We let Q =J0 
the domain Q\ 



l] 3 and choose the overlapping domain f2i t o be a rotation along the y-axis of 
= [0.3331, 0.6669] 3 as illustrated in Figure 



.1 



To examine the convergence of the 
methods, we apply the method of manufactured solutions: let u(x, y,z) = (sin(7ry) sin(7rz), 0, 0) 
and p(x, y, z) = cos(7rx) + 1. The right-hand side / is defined accordingly so that ( 1.1a I is satisfied, 



and the corresponding Dirichlct boundary conditions are strongly imposed on the entire boundary 
df2. The numerical approximation corresponding to (3.7) is then computed on a sequence of 



overlapping uniform meshes {(7~^ ,T 2 N )}% =0 - The mesh size of the initial meshes Tq and is 
ft-max ~ 1/3 and each subsequent pair , ) is generated from the previous one by uniformly 
refining each mesh. 

The stabilization parameters were set to 7 = 10 and 5 — 0.05. To solve the linear system of 
equations, we apply a transpose-free quasi-minimal residual solver with an algebraic multigrid pre- 
conditioner, filtering out the constant pressure mode in the iterative solver. The solves converged 
to a tolerance of 10~ 8 in between 34 and 47 iterations. 



The resulting errors for the sequence of refined meshes are given in Figure 8.2 Theorem 6.1 
predicts first order convergence for both the i/^-norm of the velocity error and the _L 2 -norm of 
the pressure error. This is also (at least) observed in the numerical results, thus corroborating the 
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• • A: 1.00 | 

















| • • A: 1.49 



Fig. 8.2 Convergence of the stabilized Nitsche overlapping mesh method for V h X Q^. The legend gives the fitted 
slope for each configuration. Top: Jf^-error \\u — Uh\\l,T*UT2 f° r * ne velocity versus maximal element diameter 
ft max . Bottom: L 2 -error \\p — Ph\ \j-*uT2 f° r * ne Pressure versus /i max . 



theoretical estimate. We additionally observe that the pressure approximation seems to converge 
at a higher order (« 1.5) for the range of meshes investigated here. 



8.2 Condition number tests 



The next numerical example demonstrates that the condition number of the matrix A defined 



by (7.1) is bounded independently of the position of the overlapping mesh relative to the back- 



ground mesh. This bound can be attributed to the term (3.4) defined on the overlap region: we 
also illustrate that the condition number seems unbounded if this term is not included. 

Let fl — [0, l] 3 be tessellated by iV 3 cubes, each cube divided into six tetrahedra. Define the 
overlapping domain by !?2 = ^a(0 = [1, 1 — for a parameter / E (0, 1), and tessellate Q-i in the 
same manner as f2, but with M 3 cubes. We will consider two cases of mesh sizes: (i) N = 5, M = 3, 
and (ii) N — 10, M = 6. For both cases, / approaching 0.2 from above corresponds to a limiting 
case: a few of the degrees of freedom in the overlapped mesh will only receive a contribution from 
the integrals defined over the overlap region while the cell and interface contributions vanish. Here, 
we therefore examine / S (0.2,0.21]. For each I and choice of (N,M), we compute the condition 
number k of the corresponding matrix A as the ratio of the absolute values of the largest (in 
modulus) eigenvalue and the smallest (in modulus), nonzero, eigenvalue. 

Let 5 = 0.05 and 7 = 10.0 as before. The resulting condition numbers, scaled by the square 
mesh size of the overlapping mesh, are given in Table [8T| We observe that for each I, the difference 
in the condition number between the two mesh sizes is small, as expected in view of Theorem |7.1| 
Moreover, for both the case N = 5 and N = 10, the scaled condition numbers seem clearly 
bounded as I — > 0.2. In contrast, the scaled condition number grows significantly as I — > 0.2 when 
the overlap integrals Sh are not included. 



Case 


I = 0.21 


0.201 


0.2001 


0.20001 


0.200001 


N = 5 (h ft: 0.33) 


1076 


1207 


1220 


1222 


1222 


N = 10 (h ft 0.17) 


955 


1149 


1170 


1173 


1174 


N = 5 without Sh 


583 


643 


958 


9715 


110636 



Table 8.1 Scaled condition numbers re/1 2 where h is the minimal cell diameter of the overlapping mesh. Each 
column to corresponds one overlapping domain fl(l) for I = 0.21, I = 0.201 etc. The bottom row corresponds to the 
matrix induced by the form without the overlap integrals s^. 
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8.3 Flow around an airfoil in a channel 

Finally, we illustrate that the method and technology developed here can be successfully applied 
to more realistic numerical simulations. In particular, we consider the flow around an airfoil in a 
channel. 

As 72 we take a tetrahedral mesh approximation of a sphere surrounding an airfoil, discretizing 
the boundary layer around the airfoil with a higher mesh resolution than the remainder of the 
domain. As To we take a tessellation of [— 3,3] 3 \0 where the interior domain O is contained in 
the convex hull of 75- (In total, 7o U T2 is a mesh of a [—3, 3] 3 box with an airfoil removed.) 

The stabilization parameters and choice of finite element spaces are as before: x Q\, (3 = 1, 
6 = 0.05 and 7 = 10. We enforce a parabolic velocity profile at the inflow boundary, no slip 
conditions for the velocity on the outer walls and the airfoil boundary, and stress-free boundary 
conditions at the outflow boundary and take / = 0. 

The flow patterns around the airfoil can now be studied for instance for different angles of 
attack by rotating T 2 while keeping 7o fixed. The velocity and pressure approximations for a series 
of mesh pairs, in which T2 was first rotated around the z-axis with angle 6 and then around the 
y-axis with the same angle, are visualized in Figures |8.3| and |8.4| We especially note the smooth 
transition of the solution from To to T2', the interface is not visible. 

9 Conclusions 

We have presented and analyzed a general class of stabilized finite element methods for the Stokes 
problem posed on overlapping, non-matching meshes. The interface conditions are enforced using 
Nitsche's method, yielding a provably optimally convergent method. In addition, a least-squares 
term on the overlap region is included to control the condition number of the stiffness matrix. 

The theoretical results have been verified numerically for a test problem consisting of Stokes 
flow through a channel described by a sphere superimposed on a tetrahedral mesh of the unit 
cube. We have further verified that the condition number of the stiffness matrix remains bounded, 
independently of the location of the overlapping mesh relative to the background mesh. Finally, 
we have demonstrated the applicability of the proposed method and our implementation to the 
simulation of Stokes flow around a three-dimensional airfoil with a fitted mesh superimposed on 
a non-matching background mesh. 

A natural extension of the work presented in this paper is to the simulation of fluid-structure 
interaction problems where a fluid mesh of the surroundings of an elastic body is superimposed 
on a background fluid mesh as in Figure [TTT) The fluid-structure interaction may then be handled 
via a standard arbitrary Lagrangian-Eulerian (ALE) approach on the overlapping mesh, while the 
Nitsche overlapping mesh method analyzed in this paper is used to enforce the interface conditions 
across the fluid-fluid boundary. We explore this technique further in ongoing work. 

While all software used in this work is available as free/open-source, additional work is required 
to create interfaces and documentation that make the software useful to a general audience. This 
issue will be addressed in the near future, with the goal to provide an easy-to-use programming 
environment for overlapping mesh methods as part of the FEniCS software suite. 
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Fig. 8.3 Velocity streamlines (left) and pressures in the plane denned by the z-axis (right) for different angles of 
attack Q: from top to bottom: = 20°, 0°, -20°. 



A stabilized Nitsche overlapping mesh method for the Stokes problem 



23 




Fig. 8.4 Velocity streamlines and magnitude in a cross-section for angle of attack 8 = —20°. Top: The solution 
on the background mesh T* . It can clearly be seen that the streamlines stop close to the interface zone. Bottom: 
The solution on the overlapping mesh Ti superimposed on the solution on the background mesh T* , indicating a 
smooth transition of the solution from 7~j* to 75- 
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